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INTRODUCTION  TO  VOLUME  II 


This  volume  consists  o£  a  series  of  appendices  whose  contents  could 
not  be  included  in  the  body  of  the  report  without  breaking  the  continuity. 
The  appendices  describe  topics  that  were  studied  during  the  course  of  the 
current  contract  but  are  only  tangentially  related  to  the  theme  of  Volume 
I.  Several  references  were  made  in  Volume  1  to  appendices  in  this  volume 
for  the  details  of  specific  methods. 

The  notation  introduced  in  each  appendix  applies  only  to  that  append!^ 
unless  otherwise  noted. 
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APPENDIX  A 


THE  METHOD  OF  ITERATIVE  GENERALIZED  LEAST  SQUARE 


In  Prony’s  method  we  have  the  difference  equation 
N 

y  «*  I  -  0  ,  K  -  0,1,..., Y-l  (A.l) 

pt0  P  P4* 

Y  -  M-N 


This  is  usually  solved  as  an  inhomogeneous  equation 
N-l 

ap  *"p+K  "  -IN+K  ’  K  “  0, Y-l  (A. 2) 


where  was  set  to  one. 


However  we  do  not  know  the  I^'s  exactly.  The  measurements  of  the  1^'s  are  Y. 

as 


\’h+\ 

(A. 3) 

eU 

where  eR  is  the  error  in  the  K  sample. 

Hence, 

N 

N 

J  a  Y  - 

>  a  e  .  _ 

(A.4) 

p-0  p  p+k 

pto  >  **■ 

Rewriting  the  above 

gives 

N-l 

N 

I  o  Y  - 

p-0  P  P+K 

'Y»«  +  pi0  “p  Vk  • 

(A. 5) 

N 

Letting  WR  -  ^  op 

®p+K 

(A. 6) 

A-l 


/y*^jtssi<?d£££sa&&s:-5- . 


N-l 

£  “p  W 


-Y  +  U 

V-K  K 


(A. 7) 


The  are  che  residuals  but  from  (A. €)  the  residuals  are  correlated  and  a 
least  squares  solution  to  U.6)  will  give  biased  estimates. 


When  we  use  leant  squares  to  minimize  che  residuals  by 

d-  X  wk2-0  • 

m  K^O 


ve  obtain  the  final  expression 

N~1  Y-l  Y-l 

Vci15yy«-'S'yy 

£  p  ^  P+K  m+K  £  N+K  m+R  ,  m-0 . N-i 

p»0  F  K-0  r  E>0 


We  will  get  biased  estimates  for  the  a  . 

P 

One  way  to  correct  this  is  to  use  a  method  known  as  iterative 
generalized  least  squares.  Rewrite  (A. 7)  as 


(A.  8) 


(A. 9) 


X  a_  Y  »  K"0»  1» « *  •  »Y“1  • 

pt0  P  P+K  K 


(A.  8) 


Now  define  new  notation  for  the  sake  of  convenience.  First  Introduce  the  shift 
operator  q  defined  so  that 


qf  •  f 
q  K  K+l 


q  f  m  f 

4  k  K+: 


(A.  9) 


The  polynomial  operator  A(q)  is  defined  as 


A(q)  -  f  a  q® 
n?0  ° 


(A. 10) 


Now  (A. 8)  can  be  written  as 


A(q)  Yr  -  WR  ,  let  a  -  a  .  (A. 11) 

In  matrix  form  (A. 8)  is 


Now  let  us  operate  on  WR  with  the  polynomial  filter 

B(q)  \  m  \  (A. 12) 

to  give  nR  which  is  an  uncorrelated  noise  sequence  (random  variables)  so  that 
B(q)  A(q)  YR  -  . 

Let  A  commute  with  B  so  that 

A(q)  B(q)  YR  -  nR  (A.13) 


A- 3 


and  define 


\  -  B(q)  Yx  (A. 14) 

so  that 

A(q)  Yr  -  nR  .  (A. 15) 

Now  equation  (A. 15)  has  uncorrelated  residuals  and  hence  the  solution  should 
converge  to  unbiased  estimates. 

The  iterative  procedure  is  therefore  as  follows: 

1.  Solve  equation  (A. 11)  using  the  normal  least  squares  procedure. 

This  gives  an  estimate  of  A(q)  which  can  be  thought  of  as  the 

A 

ci's  of  (A. 8).  Call  that  estimate  A(q) . 

A 

2.  Substitute  A(q)  into  (A. 11)  to  produce  an  estimate  of  the 

A 

residuals  WR. 

A  A 

3.  Use  Wg  in  (A. 12)  to  make  a  least  squares  estimate  of  B(q) . 

A 

4.  Calculate  Yr  in  (A. 14). 

A 

5.  Use  the  YR  to  make  a  new  least  squares  estimate 

6.  Continue  this  procedure  until 

2 

2,  W  no  longer  decreases  with  the  next  iteration. 

p-0  K 

The  question  now  is: 

A 

How  does  one  obtain  B(q)  of  step  3  above? 

Assume  we  have  used  the  least  squares  process  once  to  find  the  coefficients  a  . 

P 


A-4 


Next  calculate  the  y  residuals  W  as 

& 


N 

Z  °o  ^n+K  "  *  K"0*1* . . . ,Y“1 

p.Q  p  P+K  K 


The  results  of  steps  1  and  2  above  gives  a  y  diaenHional  vector  of  the  re¬ 
siduals  w  . 

Calculate  the  autovariance  function  of  the  by 
r(u)  «  ,  u-0,1 . y-1 


where 


1  Y-u 

c(u)  ■*  I  (Wi  -  W)  <W1+u  -  W)  ,  u-0,1 . Y-1 


»  -  7  I  ■ 

y  i-1 


Hence  the  autocovariance  function  r(u)  is  defined 

Y  Wt-  5;  -  5) 

r(u)  -  - -  ,  y-0,1, . . .  ,y-1 

V  (w4  -  w) 2 

i-i  1 


We  now  have  a  vector  R,  y  long. 


We  want  to  test  this  R  vector  for  whiteness  of  the  residuals  W^, 
The  standard  deviation  of  a  single  autocovariance  function  estimate  is 


The  95%  confidence  limits  for  a  single  autoccvariance  sample  are 


r(u)  +  i,96 

ST 

1  96 

Now  supposedly  if  95%  of  the  samples  r(u)  are  less  than  ■=  —  then 

Sy 

we  have  95%  confidence  that  the  noise  is  white. 

If  the  r(u)  flunks  the  whiteness  test  which  will  happen  the  first  few  times 
we  Iterate  we  must  then  go  back  and  whiten  the  residuals  Wtf. 

From  Equation  (A. 12)  we  want 

B(,)  wr  *  \  • 

L 

Remembering  that  B(q)  -  J  b  q  , 

JiV  ® 


As  an  example,  assume  L  ■  2  and  K  *  3.  Then  we  get 


W0  Wl  w2 

W1  W2  W3 

w2  w3  w4 

Wo  w.  W, 


Solving  Wb  -  n  for  the  bQ  by  using  least  squares  gives 


T  T 

WAWb  -  W  n 


b  -  (WTW)  WYn  . 


Therefore  the  new  Y_  previously  defined  as 


can  be  written  as 


h  "  Vk  +  Vk+i  +  ■'  •"  +  V: 


K+L 


APPENDIX  B 

EFFECT  OF  INCREASING  MODEL  ORDER  IN  PRONY'S  METHOD 


A  study  of  the  effect  of  increasing  model  order  in  Prony's  method 
has  been  made  and  the  results  are  presented  here  along  with  some  preliminary 
conclusions.  In  past  investigations  of  Prony's  method  it  has  beeu  noted 
that  increasing  the  model  order  above  the  known  order  of  the  waveform 
improves  Prony's  method's  ability  to  estimate  accurately  the  true  poles. 
Although  the  accuracy  of  the  pole  estimates  is  improved,  a  side  effect 
of  this  procedure  is  the  problem  of  distinguishing  between  true  poles 
and  those  poles  that  simply  fit  to  the  noise  and  have  no  relation  to  the 
Information  that  is  to  be  extracted  from  the  waveform. 

In  this  study  we  relate  the  inaccuracy  or  bias  of  the  parameter 
estimates  to  the  amplitude  of  the  residuals  of  the  least-square  Prony 
procedure.  We  assume  that  the  inhomogenous  solution  (defined  in  Volume 
I,  Section  2)  is  used  to  find  a  parameter  vector.  The  term  "residuals" 
is  Identical  to  "equation  error".  In  this  appendix,  N  is  the  number  of 
poles  modeled  and  M  is  the  number  of  samples.  When  the  residuals  are 
zero  the  least-squares  Prony's  method  either  has  reduced  to  curve-fitting 
Prony's  method  (M-2N)  or  is  processing  a  noise-free  waveform  at  the  proper 
model  order.  In  this  case,  the  pole  estimates  are  quite  accurate  and  are 
free  of  bias.  Conversely  the  bias  in  Prony's  method  is  directly  related 
to  the  magnitude  of  i  residuals. 

Tables  B.l  through  B.9  display  the  effect  of  increasing  model  order 

at  three  different  noise  levels.  o.T  is  the  standard  deviation  of  the 

N 

noise,  is  the  average  standard  deviation  of  the  residuals  over  ten 
Monte  Carlo  runs.  Each  table  shows  the  true  parameter  values  in  the  first 
column,  the  average  parameter  values  over  ten  Mone  Carlo  runs,  and  the 
variance  of  each  parameter  in  the  third  column.  For  all  nine  cases  the 
time  window  size  is  kept  constant.  Thus,  M  ■  100  and  AT  ■  0.13  seconds. 


B-l 


The  three  pole  pairs  shown  under  the  average  column  were  obtained  by 
searching  all  N  poles  to  find  the  pairs  that  were  closest  to  the  true 
poles  in  value.  In  general,  of  the  extra  poles  not  reported  here,  N-6, 
have  residues  two  orders  of  magnitude  below  those  of  the  true  poles. 

From  the  results  of  Tables  B.l  through  B.9,  the  following  observations 
can  be  made: 

1.  Prony’s  method  seems  to  guarantee  good  results  if  oR  <  c,,. 

2.  o-  <  olT  seems  to  occur  when  4N  >  M. 

That  good  results  begin  at  oR  ■  e.,,  is  not  surprising.  In  this  case  one 
would  expect  that  the  residuals  are  beginning  to  approximate  the  noise  and 
that  the  "modeled"  waveform  approaches  the  uncorrupted  waveform.  But  the 
observation  that  cR  ■  approximately  when  4n  ■  M  has  no  obvious  explan¬ 
ation.  This  observation  would  obviously  not  hold  true  if  we  were  to  apply 
Prony’s  method  to  a  noise-free  waveform.  In  this  case,  the  residuals 
would  drop  to  zero  as  soon  as  the  model  order  reached  or  exceeded  the  order 
of  the  waveform.  We  might  surmise  than  that  thir  occurrence  at  4N  -  M 
is  attributable  to  the  nature  of  the  uncorrelated  noise. 

There  is  no  simple  explanation  for  improved  accuracy  at  higher 
orders.  At  least  two  factors  seem  to  be  involved: 

1.  Increasing  the  length  of  the  parameter  vector  (by  increasing 
the  order)  is  an  effective  treatment  for  the  dense  sampling 
problem  (which  is  described  in  Volume  I,  Section  3  of  this 
report.  ) 

2.  Making  the  data  matrix  more  nearly  square  must  reduce  the 
magnitude  of  the  residuals.  Smaller  residuals  Implies  smaller 
bias  in  the  parameters. 
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Table  3.1. 


Modal  order  study  M-100,  N-12,  o  «. 
5r-.1808.  N 


TRUE  AVERAGE  a2 

-8.2000006-02  * -6.  8  129081-01 - 775463256- 03' 

>8*2000006-02  -6*  8129Q8E-01  7.546325E-03 

-1.470000E-.01  -6*  6420288-01  5.821274E-03 

-1.470000E-01  •  -6*  6420286-01 - 5T8 212746*03' 

-1*  880000E-01  -3.  3292166-01  2. 3713706-03 

-1*8800006-01  -3.3292166-01  2*3713706-03 


“97  26  0  0  00  6-01  0*  OT - 

-9*2600006-01  0*  0* 

2*  874000E+00  2*  34107664*00  3*2397366-03 

- 2. 8740 00 E ♦  00 . -2.  3416766 ♦  00 - STHWCTCTjr 

4*  8350006400  4.  84431164*00  6*0474686-04 

-4*  S3S0Q0E+  00  -4.8443116400  6*0474686-04 


'  l .  0  0  0  CO 0 6 4  00 . 1.1093956 400 - 7750'9nr36*D2' 

1*  OOOOOOE+OO  1*1093956400  7.5091936-02 

1*0000006400  1*1443786400  1*3819876-02 

- l'i  0  0  0  0 0 0  6 4  0  0 - 1.T44 3786VOO - 173  mT7 6=62' 

1*0000006400  1*1899036400  6*4747296-03 

1.000000E400  1.189903640Q  6.474729E-03 


”4*  4191116*15 - 470776926*28' 

'4.  4191116-15  4.  0  776  92E-28 

1*5079676-01  1.5089536-02 

-1*  5  079676-01 - r75TNr9'53E*Q2 

-1.  7015966-01  5*  4768  87E-03 

1.  7015966-01  5*4768876-03 


100, 


REAL  PARTS 
OF  POLES 


IMAG.  PARTS 
OF  POLES 


MAGNITUDES 
OF  RESIDUES 


RADIAN  PHASE 
OF  RESIDUES 


Table  B.2. 


Model  order  study  M-100,  N-20,  aN«.100v 
aR".l24. 


TRUE 

'T.20  0  000E-fl2 
•8.200000E-02 
■  1. 470  OQOE-Ol 
•1.470000E-01 
■1.88 OOOOE-Ol 
■1.880000E-01 


AVERAGE 

■9. 069950E-02 
•9. 369930E-02 
■1.  97178_3E-pJL_ 
•1.  971 783E-Q 1 
-1.991637E-01 
■1. 9916376-01 


6.  905491E-03 
6.905491E-05 

_ 1  •  0  BJ3JL  5,5  E^O 't 

1.0831 55E-04 
1.278673E-04 
1.278873E-04 


REAL  PARTS 
07  POLES 


UTm  00  0  E - 01  9.  2  89997 i- 01 

-9.280Q00E-01  -9. 289997E-01 

2.874000E+00  2.8832SlEt00 

-2. S7400QE  +  00  “2 .8832 
4.83500QE*00  4.8402I1E+00 

-4.833000E+00  -4. 840 21 1E*00 

■  mm  .  .  »  •  IM**)..**  I  .  “•  -  «>«•«■>•  m  ■  ••  n  m n  —  ■ . ie^—i 


3.6S7123E-05 
3.657123E-05 
1*  0110442-04 

TTiTfnOT=Tf 
1. 1760392-04 
1^1780392-04 


IMAG.  PARTS 
OF  POLES 


l.  000  0 o' 0 E ♦' 0 0  '  “ ' 7 U  027 8 5 9 E ♦  00  2.  0 3744 3E- 03 
1. OOQOOOE  +  OO  1.  0  27 83 9E*0 0  2.037443E-03 

...l*  M.WQ0  E  top. _ 1.0.?_84&71±jlQ _ fcfli 

1. 0  0QQQ0E  +  00  1.  0364S7E+00  2.199189E-03 
l.OQOOOOEfQO  1«  018026E+00  1.059788E-03 
1 . 0 0 0 0 0 0 E ♦  00 if  018.0 ? 8 g*D0  1.  Q997.88E-03 


MAGNITUDES 
OF  RESIDUES 


0. 

-1. 230493E-02 

9.820293E-04 

0. 

1. 230493E-02 

9. 820293E-04 

0. 

-1. 9287782-02 

2*  0420 11E-03 

0. 

lelZSftfC^OZ 

2.0420 11E-03 

0. 

-9*  326322E-03 

2.024303E-03 

0. 

9.  328322E-03 

2.  0  243  03E-03 

RADIAN  PHASE 
OF  RESIDUES 


B-4 


Table  B.3.  Model  order  etudy  M-100,  N-36,  o„-.100, 
3„-.0678. 


TRUE 

8.200000E-02 
8*  2000002*02 
1.470000E-01 
1.470000E-01 
1.  8800001*01 
1.880000E-01 


AVERAGE 

ei  1691  ME- 02 
•6.169179E-02 
■1.  483J3QE-Q1 
•1.4839381-01 
•1.  842S90E-01 


3. Q56697E-03 
3.096697E-05 
.V07.22JJE-.QA 
!• 1722  73E-04 
1.861397E-04 


REAL  PARTS 
OF  POLES 


9.  26 0  0 OOf- OT 
9.  260000E-01 
2.S74000E+00 
Z.874000E+00 
4.  639000E+  00 
4. 839000E+00 


•9.26648  QI-01  2.394936E-Q9 

2.8748361+00  9.301869E-05 

•  2.8  74  8361  ♦ OH - 9731)1 8'6TE-0J' 

4. 836128E+Q0  4.490036E-09 

•4.8361281+00  4.490036E-09 


IMAG.  PARTS 
07  POLES 


1. OOOOOOE+OO 
1. 0  OOOOOE  +  00 
L.OOOOOOE+OQ 
L.  OOOOOOE  +  OO 
1. OOOOOOE+OO 
1. OOOOOOE+OQ 


9.931078E-01  1.329078E-03 

9. 93107BE-01  1.329078E-03 

..1*  0137201+00  4. 119833E-03 

1. 01372  bE  +  00  4. 119833E-03 
9. 691063E-01  6.276799E-03 

9.  6 91 0 » 3E-01 6x2767991-03 


MAGNITUDES 
07  RESIDUES 


•3.’963147E-oT 
3. 963 147E-03 
1.  449398E-03 
•1.44639  8  Ml 
•1.  402673E-02 
1.  402673E-02 


879139991-04 
8. 993999E-04 

iiLl9UlbA.l 
177902291-03 
7. 1 34940 E- 04 
7.134940E-04 


RADIAN  PHASE 
07  RESIDUES 


Tfebla  B.4.  Modal  ordar  itudy  M-100,  N-8,  c^-. 00100, 
oR-. 007098. 


TRUE 

AVERAGE 

a 

-TT2TT  tnjm-TTi 
-8.2000008-02 
-1. 4700008*01 

-5. 7480i Jg-oi 
-4»  746 0216-01 
-3.23707&E-01 

^.f7l5TO?fl3 
6.  8715  08E-03 
1.1543B3E-03 

•1.  47Q0tttJE-()i 
-l.seooooE'Oi 

-1.880000E-01 

-yrrynnrnnr” 

-1.9949996-01 
-1. 9949338-01 

“TTmuTJi^ir 

2. 6146S7E-P5 
2. 614657E-09 

• 

9«  2600006-01 
-9.260000E-01 
2* 8740008400 

8.9503058-01 
-8. 9503881-01 
2.8743348400 

8.7045 12E-04 
8. 7045 126-04 
2.6016698-04 

-i.iHoooIioo 

4.8350QQE+00 
-4.835 OOOEfOO 

-2.8^43941400 

4. 843241(400 
-4. 843241E400 

2. 60lZ  B5E-0  + 
3.f 698681-05 
3. 669866E-0! 

I.OOOOOOE+OO 
1.  OOOOOOEMJO 
1.  000  QQOE4-00 

177371 8TE4  00 
1.7371891+00 

1.  223620E+00 

9.1143  S6E-0T 
9.1143668-03 
1.4617208-03 

l.  doooobEiofl 
l.OOOOOOE+OO 
1.  00 0 OOOE <*•  00 

ummrm— 

9. 0283708-01 

9.  02837  0e-01 

1.461720E-0S 
4. 5500  50E-04 
4.5500508-04 

REAL  PARTS 
OF  POLES 


IMAG.  PARTS 
OF  POLES 


MAGNITUDES 
OF  RESIDUES 


"0.- 

2.0228131-01 

“77ZT55  541-03 

0. 

-2.0228038-01 

7.2355S46-03 

0. 

2.2963318-01 

3. 7694  24E-03 

0. 

-2. 29B3llX-01 

3.7694X4E-03 

0. 

1.2608948-01 

2.4557  851-05 

0. 

-1.2806948-01 

2. 455765E-0S 

RADIAN  PHASE 
OF  RESIDUES 


Table  B.5.  Model  order  study  M-1G0,  N-12,  a 


TRUE 

3R-. 00300. 

AVERAGE 

A 

-8.2000001-02 

-6.3347901-02 

3.6196361-08 

-6.2000001-02 

•6. 3347901-02 

3.  6196361-06 

-.1  ..4,7-0  0  0  0  EiO  1 _ - 1  ..4620  4i  £s.01_. 

_ 3  ..10  7  9  9  6 1.?.0  6. 

-1.4700001-01 

-1.4820661-01 

3.307596E-08 

•1. 66  OOOOE-Ol 

-1.679962E-01 

1. 695690E-07 

| 

I 

• 

_ lju63.949.fl.l-JI. 

9.260000E-01  9.  2632271*01 

-9.2600001-01  -9.2632271-01 

(*.»  839 000 E4-  00  4.  6390411+00 

-4.6390001+00  -4. 6390411+00 


3. 0942061-06 
3. 0  942  Q6E-06 
2.3133361-06 

iltimaHi 

6.2970911-02 

8.7970911-02 


1.0000001+00  1. 0094411+00  9.2241061-09 

1. 00  QQQQE  +  00  1.  0094411+00  9.2241061-0$ 

_l.-OJUU.Ofl It 00  .  .l..Jia3AgJltJfl  —  4.7-160  98 t-JJL 

1.0000001+00  1.0036931+00  4.2160961-09 

1. 0000001+00  9.9943121-01  4.0640301-06 

1. 0000001+00  ,9.  994  3121-01  4.  0640  301-0  Bl 


0* 

-1.1104161-03 

6.0193231-09 

0. 

1. 1104161-03 

6.0193231-09 

0. 

6.  180  273 1-04 

4.0663171-09 

0. 

-8.1802701-04 

4.0663171-09 

0. 

1.  6991911-03 

9.1137111-06 

_0a _ 

-i.  »mi[-3i 

9.1137111-06, 
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-.00100, 


REAL  PARTS 
OF  POLES 


IMAC.  PARTS 
OF  POLES 


MAGNITUDES 
OF  RESIDUES 


RADIAN  PHASE 
OF  RESIDUES 


Table  B.6.  Model  order  study  M-1G0,  N-20,  aN-.001, 
aR-. 001246. 


TRUE 

•8.200000E-02 
*8*20  OOOOE-02 


AVERAGE 

'6«  199835E-Q2 
•8.199885E-02 


6.959435E-09 

6.959435E-03 


•  A|  f  f  W  »  ■  ^  •  W  * 

-1.470000E-01 

-1.580000E-01 

-1.580Q00E-Q1 

+  .T.  ww#.,  y  * 

-1.  47013SE-Q1 
-1.8S03S8E-Q1 
>1.8803982-01 

1*  0  34647S-08 
1. 121223E-08 

9.260000E-01 
-9. 260000E-01 
2> 8740QQE+00 

9.  2600ME-01 
-9.26007BE-01 

2. 874021 EfO  0 

3. 504264E-09 
3.  504264E-09 
9. S360S2E-09 

•2.874000E*00 
4* 835000E+  00 
-4._83EQ_0Q_L»M_ 

-2*  S74025E+00 

4.  834994E+00 
-4.834997E+Q0 

9.536052E-03 
1.374833E-08 
1. 374833E-Q8 

1.  000000EM2 
it  OQOOOOE+OO 
l.OOOOOOEfOO 

9. 999503E-01 

9.  9999 B9E-Q1 

1*  OQOllOE  +  OO _ 

2.130885E-07 
2. 1308  8SE-Q7 
..  2i_0.91818E-0r 

t.OOOOOOEfflO 
1*  QOOOOOE>QO 
1.000000E+00 

1*  OOOllOE+OQ 

1.  000 088E  +  0  0 

1.  000088E+00 

2. Q91618E-Q7 
9.094737E-08 
9. 094737E-08 

inr 

0s 

0. 

-47  4  Q150ZE-0E 

4. 401502E-0S 
-8.  97328  0E-0E 

9. 03S431E-T8 
9.035431E-08 
2. 0  04572E-Or 

0. 

0. 

JU _ 

8»  9732SOE-0S 
-2#  903918E-0E 
2.  903918  E-Q9  _ 

2.  004972E-Q7 
2. 259596E-Q7 
■■JiillHSAE zAL 

REAL  PARTS 
OP  POLES 


XMAG.  PARTS 
OF  POLES 


MAGNITUDES 
OF  RESIDUES 


RADIAN  PHASE 
OF  RESIDUES 
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Table  B.7.  Model  order  study  M-100,  N«8,  a.,-. 0100, 
or-.0394. 


TRUE  AVERAGE  a2 

8.200000S-Q2  -4. 81S941E-01  4.313874E-Q3 

8*  20  Q0Q0E-02  -4* 81694LE-Q1  4.313874E-03 

1*470000E-01  *4i 931 0196*01  1.633214E-03 

1.470000E-01  -4.931019E-Q1  l*6332UE-03 

1.88000QE-Q1  -2.603548E-01  4. 5 l 87 06E-05 


REAL  PARTS 
OF  POLES 


9.260000E-01 
9.260Q00E-01 
2.87400QE»0Q 
2.874000E4QQ 
4. 835000E<*  Q0 


■2.444876E+00 

4.820967E+00 


6.252300E-04 
1*6992  93E-04 


IMAG.  PARTS 
OF  POLES 


1*000000£»00  8.254609E-01  2.9306S2E-02 

l.OOOOOOEtOO  8.234609E-01  2*9306821-02 

1*OOOOOOE»OQ  1* 194068E+00  3*368307E-03 

1. 00 0000 S ♦bo  l.tiVoV8E400  3*i68307E-Q3 

1*  OOOOOOE  +  QO  li  l 12501 E+00  4.643365E-05 

1*  0000  g OE  » 0 0  1.112301E+00  4, 643369E-05 


MAGNITUDES 
OF  RESIDUES 


2*  991933E-15 
1*  623  083  E-01 
•1.6250S3E-01 
•l*  073143  E-01 


6*  5  96341E-28 
2.  7298  88E-03 
2*  729888E-03 
9*  082337E-04 


RADIAN  PHASE 
OF  RESIDUES 


'l 


I 


I 


ia 


Tible  B.8.  Model  order  study  M«100,  N**12,  g  «.0100, 
or-.0295. 


TRUE  AVERAGE  g2 

8 . 2  0  0  0  0  0  E  -  02  -"2*  353161E-01  2  a  3 5  92  6  8E- 0 3 
8*  200QQQE-02  -2. 35315 IE- 0 1  2.859268E-03 

1*  470  POPE*  01  -2.  4074SS£*01  5.513517E-04 

1.470000E-01  -2. 407455E-01  5. 5 1 S517E-04 

1.880000E-01  -1.934581S-01  1.094043E-05 

1 . SHOOOOE-Ol  -1. 934551E-01  l«094Q43E-05 


REAL  PARTS 
OF  POLES 


9*  26 0 Q0 0S* 01 
2*  87400QE+  00 

9.4d447J£-di 
-9*  4 0887  9E-01 
2.  899695  E+  00 

3.593158E-04 

1.687968E-04 

r2TST4  000  ETO  fl~ 
4*  8350002*00 
>4*  8350002  +  00. 

"^7T9f9¥37ETTir 
4*  842363 E+00 
-4*  S423S3E+00 

r." 617  9641- O'* 
1*  1226 18E-04 
1*1226182-04 

1.000000E+00 
1* 0000002+00 
1*  OOOOOOE+OQ 

1.  4  J606QE+00 
1*  426 050 E+00 
1.2 11 07  5  E+00 

1*70  80  38E-02 
v  1*70 80 382-02 
v*  2* 967772E-03 

1*  OOOOOOE+OQ 
1*  OOOOOOE+OO 
1*  00000024-00 

1*  211 07SE+00 
9*  755319E-01 
9.755319E-01 

2*  967772E-03 
4* 1506  53E-04 
4* 1506  53E-04 

s.  Hmye-ai — r.sr^ttt-as 

2.712845E-03  7.626122E-Q3 

7»  308095E-Q2  3.766332E-03 

y^TWTOTFBZ - 3;7tT332E-fl3 

1.025553E-01  8.128753E-04 

1. 0255SPE-01  8*1287 53E-Q4 


IMAG.  PARTS 
OF  POLES 


MAGNITUDE 
OF  RESIDUES 


RADIAN  PHASE 
OF  RESIDUES 


1 


i 

I 
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APPENDIX  C 


E  ADAP 


HSilEag 


The  adaptive  method  reeulta  from  the  generalized  schema  of  Volume  I* 
Section  2,  under  the  following  assumptions: 

1.  a  “  1  and  g  *0. 

o  o 


2.  P. 


,  1  •  1, . . .n  and  F 


3.  The  model  input  is  a  unit  sample  at  k  -  0  (discrete  Impulse). 

The  unique  feature  of  the  method  is  that  the  filter  poles,  z^,  may  be 
adjusted  to  any  value  In  the  Z-plane,  An  adaptive  technique  for  adjusting 
the  filters  consists  of  first  initializing  the  filters  to  arbitrary  values 
in  the  Z-plane  and  repeating  the  following  steps: 

1.  Find  an  estimate  of  the  process  transfer  function  using 
the  current  filters  in  the  model. 

2.  Set  each  filter  pole  to  one  pole  of  the  estimated 
transfer  function. 

This  procedure  is  repeated  until  approach  zero.  The  poles  of  the 
process  can  be  estimated  during  the  course  of  iteration  by 

A  *i 

zi  1  +  oi 

removing  the  need  for  finding  the  roots  of  a  polynomial.  The  filter  poles 

are  updated  to  on  each  iteration.  When  the  approach  zero  the  pole 

updating  ceases  and  the  method  converges.  At  each  iteration,  the  a.  are 

1 

found  using  any  of  the  techniques  for  finding  a  parameter  vector  described 
in  Volume  I,  Section  2.  Perhaps  the  simplest  method  is  to  choose  the 
parameter  vector  as  the  weakest  eigenvector  of  ft*Q.  At  convert  jnce  the 
S-plane  poles,  s^  can  be  obtained  from  the  filter  poles  by 


and  the  A. 


The  adaptive  filtering  method  seems  to  be  a  very  useful  technique 
because  it  uses  what  seems  to  be  optimal  filters,  band  pass  filters,  for 
estimating  pole  locations.  In  addition  it  presents  a  measure  of  the  error 
and  iteratively  adapts  the  filters  until  the  error  is  at  the  level  where 
it  is  desired. 

As  an  example  of  the  use  of  the  adaptive  filtering  method  consider 
the  data  shown  in  Figure  C.l.  These  data  were  numerically  generated  by 
using  the  time  domain  computer  code  TWTD  [C.l],  The  structure  modeled 
by  TwTD  was  a  thin  cylindrical  scatterer.  The  signal  was  contaminated 
with  noise  to  give  a  15  dB  signal-to-noise  ratio.  Figure  C.2  shows  the 
resulting  poles  from  five  Monte  Carlo  trials. 

The  following  statements  can  be  made  about  the  adaptive  method  after 
studying  it. 

The  adaptive  method  is  a  new  method  which,  in  many  cases,  pro-  ‘des 
excellent  pole  estimates  under  difficult  conditions.  The  method  is 
unique  in  that  a  solution  to  a  polynomial  is  not  required  to  find  estimates 
of  the  process  p.les.  the  method,  in  effect,  "swallows"  the  polynomial 
solver  in  its  own  iterative  pole-searching  scheme. 

Attempts  to  analyze  waveforms  consisting  of  highly  damped  exponential 
components,  such  as  the  transient  responses  of  a  sphere,  have  not  been 
successful.  The  adaptive  method  does  not  converge  for  waveforms  which 
display  double  pole  characteristics,  that  is,  waveforms  with  components 
of  the  form  t  exp(st).  Slight  modifications  to  the  adaptive  method  might 
allow  the  analysis  of  such  waveforms. 
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APPENDIX  D 


THE  PENCIL-OF-FUNCTIONS  METHOD 

The  pencil-of- function*  method  results  from  the  generalized  modal 
daacrlbad  in  Voluma  I,  Saction  2  with  Cha  following  aaaumptiona: 

1.  Oq  ■  1  and  8q  ■  0. 

2.  (a )  ■  (1/a)*  (caacadad  concinuoua  intagratora) . 

This  mathod  ia  not  vary  aaay  to  implement  on  a  digital  computar  ainca  tha 
continuous- time  intagratora  cannot  ba  imp lamantad  exactly  by  any  algorithm. 
Whan  approximate  intagratora  are  caacadad,  aa  thay  ara  in  the  pancil-of- 
functiona  mathod,  largo  arrora  can  bo  quickly  accomulatad  and  tha  intended 
raault  after  a  number  of  intagrationa  daatroyad. 

Thia  difficulty  can  ba  raaolvad  by  caacading  diacrata  Intagratora  to 
obtain  filtara  with  pulaa  tr ana  far  functiona  given  by: 

rt  *  'i(,>  ■  (A)  ‘  (if 

which  can  ba  lmplemantad  on  a  digital  computar  with  no  error  by  ualng 
difference  equatlona.  Tha  variable  Z  la  dafinod  aa 


2  "  1  a 


and  z  is  tha  z-tranaform  variable. 

Whan  diacrata  intagratora  ara  uaad  the  diacrata  pancil-of-functiona 
mathod  raaulta.  Tha  polaa  of  tha  pulse  transfer  function  of  the  process 
or  wavaform  mey  ba  estimated  aa 


si  “  i-V 


kL 

where  is  ths  i  zaro  of 

zn  +  a-  Zn-1  +  ...  +  an  -  Q  (D.l) 

l  n 

In  the  original  pancil-of-functions  mathod  tha  vara  found  by  using 
Jain' a  mathod  for  constructing  tha  parameter  vactor  which  la  daacrlbad  In 
Voluma  I,  Sactlon  2.  With  Jain's  mathod 


aU  ak 

whara  A..,  in  this  case,  is  tha  olamant  at  tha  1  row  and  j  column  of 
•dj  fl*n.  Othar  mathods  can  ba  usad  to  astimata  tha  parameter  vactor. 

Tha  S-plane  astimatas  of  tha  poles,  s^,  ara  ralatad  to  tha  zaroa  of 
(D.l)  by 

lnO-Z£, 

I 

Ona  difficulty  with  tha  pancil-of-f unctions  mathod  Is  ralatad  to 
tha  ar.tanuation  of  tha  highar  fraquancy  modas  of  tha  procaaa  output  by 
tha  rapaatad  integrations  applied  to  tha  output  waveform.  It  can  ba 
verified  that  an  intagrator  is  simply  a  first  ordar  filtar  whosa  Laplaca 
transfar  function  has  a  polo  at  tha  origin  in  tha  S-plano.  Such  a  filtar 
tands  to  suppraas  tha  highar  fraquoncias  prasant  at  its  input.  Tha 
highar  fraquancy  supprassion  phanomanon  is  lllustratad  in  Figura  D.l. 
Normally,  whan  an  axponantial  function  is  intagratad  rapaatadly, 
componants  of  powar  of  tlma  axlst  in  tha  highar  intagrals  as  wall  as  tha 
original  axponantial  function  componants.  In  Figura  D.4  tha  componants  of 
powars  of  tima  havo  baan  subtrsetod  from  tha  intagratad  wavaforms  in 
ordar  to  maka  tha  attonuation  of  tha  highar  modaa  more  avldont.  Tha 
first  wavaform  is  a  hypothatical  wavaform  provided  for  analysis.  Tha 
wavaforms  that  follow  ara  tha  intagrals  of  incraaslng  ordar  of  tha  first 
waveform  end  display  tha  incraaslng  domlnanca  of  ths  fundamantal  mode  or 
mode  of  lowest  fraquancy.  Further  integrations  yield  nearly  identical 
wavaforms.  The  intagratad  wavaforms  tand  to  become  linearly  dependant 
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wavaforma.  The  integrated  wavaforma  tend  to  become  linearly  dependant  at 
higher  nodal  ordara.  The  matrix  than  tanda  to  alngularlty  and  the 
method  becomea  unatabla.  The  auppraaalon  phenomenon  occura  in  both  the 
dlacrata  and  continuoua  mathoda.  In  fact,  even  for  vary  modaat  modal 
ordara,  the  method  can  become  numerically  ill-conditioned  to  a  degree  that 
apecial  care  muat  be  taken  to  aaaura  accurate  Invaralon  of  n*n. 


APPENDIX  E 


METHOD  OF  REDUNDANT  AVERAGING 


Tha  radundant  averaging  schema  la  a  praprocaaaing  schema  Chat  attempts 
to  combina  redundant  data  (that  la,  mora  data  than  ara  necessary  to  datanalna 
tha  parana tara)  In  a  way  that  avolda  tha  blaa  lntroducad  by  ualng  a  laast- 
aquaraa  achama  to  combina  radundant  data.  Tha  mathod  attampta  to  tranaform 
a  raw  waveform  of  aora  than  2N  aamplea  whara  N  la  tha  modal  order  Into  a 
praprocaaaad  waveform  of  exactly  2N  aamplea  by  averaging  within  tha  waveform. 
Tha  averaging  can  be  dona  ao  that  tha  axpactatlona  of  tha  polaa  of  tha  pra- 
procaaaad  waveform  ara  aqulvalant  to  tha  axpactatlona  of  tha  polaa  of  tha 
raw  wavaform  provided  tha  additive  noiaa  on  tha  raw  waveform  la  aaro  mean 
and  uncorralatad  batwaan  aucaaalva  aamplea.  Tha  praprocaaaad  wavaform  may 
than  be  procaaaad  with  cuiva-flttlng  Prony'a  mathod  to  avoid  tha  blaa  of  tha 
least-aquaraa  procedure. 

Tha  moat  general  daaorlptlon  of  tha  radundant  averaging  achama  can 
ba  atatad  almply  aa 

v1 

'  *  jlo  Vm,  ’ k"0’1 . 21,-1 

aU 

whara  N  la  tha  modal  order,  X.  danotaa  tha  k  sample  of  tha  praprocaaaad 

waveform,  and  y^  danotaa  tha  k  aampla  of  the  raw  wavaform.  In  order  to 

limit  thla  daacrlption  to  tha  aaaantlal  faaturaa  of  tha  radundant  averaging 

achama,  N. ,  tha  number  of  avaragea,  an''  N  ,  tha  decimation  epoch,  ara  not 
A  ■ 

explicitly  defined  hare.  Theae  parameters  ara  chooaan  aa  daalrad  but  usually 
in  a  way  that  would  produce  a  dealrable  praprocaaaing  mode  in  soma  sense. 

For  inatanca,  the  value  for  tha  decimation  epoch  might  ba  chosen  on  the 
basis  of  tha  maximum  frequency,  <u  ,  present  in  tha  raw  wavaform  (if  this 
information  is  available)  and  tha  number  of  averages  might  be  chosen  so 
that  every  sample  in  tha  raw  wavaform  is  used  once  tha  value  of  tha 
decimation  epoch  is  sat.  Hence  tha  sampling  interval,  At,  for  tha 
praprocassad  wavaform  could  ba  obtained  as 


i; 


4t  ■  “r« 


where  N  -  Integer  Part 


fr-ft-l- 

L  m  rawJ 


Tha  number  of  avaragaa,  N^,  can  ba  computed  by 


\  "  M  “  Na  (2N-1) 


whara  H  la  tha  numbar  of  aamplaa  in  tha  raw  waveform  and  N  la  tha  ordar 
daairad.  It  ahould  ba  notad  that  tha  abova  mathod  for  determining  and 
N(  la  not  tha  only  poaaibla  tachniqua. 


Tha  radundant-avaraging  procadura  la  aff actively  two  oparatlonas 


1.  Apply  a  low-paaa  noving-avaraga  filter,  A(z),  to  tha  raw  waveform. 

2.  Decimate  tha  filtered  wavaform  with  decimation  epoch  N  . 


Tha  tranafar  function  of  tha  moving-average  filter  las 

V1  N.-2  A  , 

A  ( m )  ■  a  *  +  a  A  +  ...  +  a+1  • 


The  numerator  of  A(r)  can  bo  factored  into 


na 


II  (Z-Z.), 
i-1  1 

The  zeroe,  z^,  of  the  numerator  of  A(z)  have  unit  magnlcudoa  and  argumanta, 

...  -  .  »»(!-» 

“*  *1  N  * 

A 

for  i  ■  1,  N^.  It  ia  clear  that  tha  firat  factor  cancela  with  tha 

denominator  giving 

Na 

A(z)  ■  H  (z-z.), 
i-2  1 
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ZEROS  OF  A ( Z ) 
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Figure  E-l.  Z-plane  plot  of 
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Table  E.l  True  poles  and  results  of  five  Monte  Carlo  runs  for  the 
redundant  averaging  example  shown  in  Figure  E.l. 


REAL  PARTS  OF  POLES 


TRUE 

POLES 


MONTE  CARLO  RUNS 


2  ' 

3 

-.127 

-.122 

-.157 

-.239 

-.103 

-.227 

AVER. 


-.107 

-.187 

-.271 


X  DEV. 


IMAGINARY  PARTS  OF  POLES 


MONTE  CARLO  RUNS 


TRUE 

POLES 


.926 

2.874 

4.835 


2 


.904 

2.913 

4.899 


5 


.936 

2.882 

4.803 


AVER. 


.925 

2.902 

4.899 


Z  DEV. 


A(z)  chen  has  zeros  which  are  evenly  spaced  around  the  unit  circle  but 
Che  zero  at  z»l  is  canceled.  Since  there  is  never  a  zero  at  Z"l,  the 
moving-average  filter  can  be  viewed  as  a  type  of  low-pass  filter,  A 
z-plane  plot  of  the  zeros  is  shown  in  Figure  E.l. 

As  an  example  of  the  application  of  this  method  consider  the  noise- 
free  waveform  shown  in  Figure  E.2a.  This  waveform  consists  of  400  data 
points  representing  a  sixth  order  exponential  function  generated  with  the 
poles  listed  in  Table  E.l.  This  data  was  then  corrupted  by  adding  whits 
noise  with  a  standard  deviation,  o,  of  0.5.  Figure  E.2b  shows  the  noise 
contaminated  signal.  This  signal  has  a  signal  to  noise  ratio  of  15.6  dB 
where  the  signal  to  noise  ratio  is  defined  as 

S/N  -  20  log 

where  Rpeail  is  the  peak  amplitude  of  the  transient  signal,  Five  Monte 
Carlo  trials  were  run  on  this  data  using  the  redundant  averaging  method. 

The  model  order  was  selected  to  be  8  (two  more  than  the  known  order)  and 
Ng  and  were  chosen  to  be  16  and  160  respectively.  Figure  3.2b  shows 
the  reconstructed  waveform  obtained  from  one  of  the  Monte  Carlo  trials. 

Note  how  good  the  fit  is  considering  the  signal  to  noise  ratio  was  15.6  dB. 
Table  E.l  lists  the  results  of  the  five  trials  and  shows  the  pole  averages, 
standard  deviations  and  per  cent  deviations.  It  should  also  be  noted  Chat 
the  signal  to  noise  ratio  compared  to  each  pole  residue  Is  only  6  dB. 

Hence  from  poles  with  a  6  dB  information  content  we  were  able  to  recover 
them  very  accurately. 

The  redundant  averaging  scheme  has  two  difficulties  or  limitations 
that  can  cause  the  method  to  be  lass  effective: 

1.  If  the  sampling  rate  in  the  preprocessed  waveform  is 
sufficiently  low,  the  higher  frequency  poles  can  be 
folded,  perhaps  several  times,  about  the  Nyquist 
frequency. 
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2.  The  method  can  null  modes  in  the  preprocessed  waveform 
that  exists  in  the  raw  waveform. 

The  frequency  folding  phenomenon  introduces  an  ambiguity  in  that  it 
is  not  known  how  many  times  a  particular  pole  has  been  folded  about  the 
Nyquist  frequency.  Hence,  N  possible  poles  are  introduced  for  each 
extracted  pole  by  the  redundant  averaging  scheme  where 

N  » 

3 


(AT)  preprocessed 
(AT)  raw 


Of  course,  if  the  highest  frequency  mode  of  the  waveform  is  known  to  be 
lower  in  frequency  than  the  preprocessed  waveform's  Nyquist  frequency,  the 
ambiguity  is  resolved  but,  in  general,  this  will  not  be  the  case. 

The  mode  nulling  can  occur  if  the  averaging  parameters  are  such  that 
the  zeros  of  the  resulting  low  pass  filter  cancel  poles  of  the  signal 
Itself.  That  is,  if  we  define  the  preprocessed  waveform  P(z)  as 

M 

P(z)  -  A(z)  D(z) 

resulting  from  the  averaging  process  A(z)  operating  on  the  signal  D(z), 
then  if  A(z)  and  D(z)  each  have  terms  in  common  they  can  tend  to  cancel 
each  other. 
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COLUMN  PRONY'S  METHOD 


i 


* 


la  this  method  the  2-plane  estimates  of  the  poles  are  the  roots  of 
the  polynomial 

^  ,  N.  ,  ,  N.2  .  t  N.N  . 

o'.q  +  0^(2  )  +  a2(z  )  +  ...  +  djj(z  )  -  0 

The  are  determined  from  the  system, 


y0  yN  *•*  yN2-N 

l 

*< 

z 

to 

_ 1 

yl  yN+l  *  *  *  YN2-N+1 

°i 

yN2+l 

•  •  • 

• 

• 

kS  • 

a 

•  •  • 

• 

• 

• 

yN-l  y2N-l  ’  *  *  YN2-1 

_°N-t 

„yN2+N-i. 

> 

where  ct^  is  assumed  to  be  one  and  {y^,  y^,  . ..,  yjj2+jj_^^  denotes  the 
sequence  for  the  waveform  containing  N^+N  samples.  It  should  also  be 
noted  that  the  polynomial  has  N^  roots  only  N  of  which  are  related  to 
estimates  of  the  true  constituent  poles  of  the  waveform. 

2 

Column  Prony's  method  offers  a  means  of  combining  N  +N  data  points 

compared  to  the  2N  data  points  of  the  standard  Prony's  method.  The 

column  Prony's  method  can  either  be  used  in  the  least-squares  version  or 

the  curve-fitting  version.  If  the  curve-fitting  version  of  column  Prony's 

method  is  used,  the  resulting  parameter  estimates  are  unbiased. 

2 

Unfortunately  the  method  yields  N  pole  estimates  only  N  of  which  are  the 
true  poles.  It  appears  then  that  column  Prony's  method  leads  to  the  same 
problem  that  Increasing  the  model  order  led  to  in  the  standard  Prony's 
method:  an  ambiguity  in  the  identity  of  the  true  poles.  Moreover,  this 
method  has  an  additional  problem:  the  matrix  can  become  nearly  singular 

I 


F-l 


( 


even  if  the  model  order  is  lower  than  or  equal  to  true  order.  This 
phenomenon  is  similar  to  the  singularity  of  the  standard  Prony's  method 
when  the  model  order  equals  the  waveform  order  and  the  highest  frequency 
is  equal  or  nearly  equal  to  the  Nyquist  frequency. 
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EVAN'S  AND  FISCHL'S  METHOD 


I 


» 


► 


» 


In  Evan's  and  Fischl'a  method  [G.l]  they  define  a  "true  error" 
sequence  {e^*  k*0,  . M-l}  as  the  error  between  the  given  waveform  and 
the  "fitted"  waveform.  They  then  proceed  to  define  the  "equation  error" 

sequence  {d^,  k*0 . M-l}.  They  further  proceed  to  define  a  relation 

between  the  equation  error  and  the  true  error: 

a  ■  Wd 


where 


•■[•q  *1  •••  Vl1  ’ 

i  3  [d0  dl  dM-N]T  ' 


and  W  ■  A[ATA]  ^ 


where 


r*  a. 


“N  aN-l 


1-0 


...  0  “1 


•  •  e  • 


0 

a0 


■  •  #  • 


aN-l 
°N  - 


M  is  the  number  of  samples,  N  is  the  number  of  poles,  0^*1,  and  the  a^; 
are  the  coefficients  of  Prony's  difference  equation  defined  in  Section  2, 
Volume  I. 


By  ueing  this  relation  between  the  two  errors  they  describe  two  iterative 
procedures  aimed  at  minimizing  tha  "true"  error  criterion: 


Ic  can  ba  shown  chat 


W(a)  - 


3a 


N-k 


A(a) 


+  A(a) 


-  W(a) 


3a 


3a 


N-k 


Av,  .) 


A(a) 


N-k 


A(a) 


A(a)4  A(a) 


and  [(3/3aN-lt)  A(a)]  is  simply  tha  matrix  A(a)  with  onaa  raplacing  tha  aN_^ 
with  all  othar  alamanes  taro. 

Tha  following  observations  can  ba  mada  about  tha  method t 


1.  It  produces  optimal  pole  estimates  in  tha  sense  that  it 
minimizes  "true  error".  This  means  that  tha  mean-square  error  between  tha 
given  wavaform  and  tha  fitted  waveform  is  minimized.  These  optimal  pole 
estimates  are  obtained  only  in  the  second  iteration  phase. 

T 

2.  .ue  matrix  A  A  must  be  inverted  on  each  iteration  of  both  tha 

T 

first  and  second  procedures.  When  the  wave .  .  ,u  has  over  100  samples,  A  A 
becomes  very  large,  and  hence,  expansive  tu  invert.  Consequently,  the 
method  is  prohibitively  expensive  when  the  wavaform  has  over  100  samples. 
This  method  has  yet  to  be  evaluated  by  tests  on  noisy  data.  Nothing  is 
presently  known  about  its  convergence  characteristics  of  its  tolerance  of 
noise. 


The  following  example  illustrates  the  optimal  estimates. 
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Exaaola 


N  •  1,  M  -  3 

2"  [2  1  2]  (vavaforo  eo  ba  approxiaatad ) 
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Initial  paramatarai  ,  a  ■-! ,  •  Thaaa  ara  tha  optimal  paramatara, 

tharafora  tha  toathod  should  convarga  ionadlataly. 
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Tharafora  tha  updatad  valua  of  la  aquivalant  to  lea  old  valua,  hanew  tha 
maehod  haa  convargad  at  tha  optimal  paramatara  for  tha  aacond  ltaration 

phaaa. 


G-5 


.... - Ua  iikt  h  il.«i J  i ki^-.  tgiuwi.-.-  ait'iidiilkj  J 


£SL£. 


Reference 


C.l  A.C.  Evens  and  R.  Pischl,  "Optimal  I *ast-Squaraa  Tims-Domain 
Synthesis  of  Racursiva  Digital  Filtara",  IEEE  Trans.  Audio  and 
Elactro-acouatiea.  AU-21,  February  1973,  pp.  61-65. 


06 


APPENDIX  H 

CONSTRAINED  PRONY 'a  METHOD 
H.l  Uni  of  a  Constrained  Prony  Method 

Many  tlmea,  In  tha  application  of  tha  Prony  algorithm  aoma  of  tha 
polaa  ara  known  a  priori.  It  would  ba  vary  uaaful  if  tha  Prony  algorithm 
could  ba  conatralnad  in  auch  a  mannar  that  tha  knowladga  of  tha  known 
polaa  la  uaad  in  axtractlng  tha  unknown  polaa.  Tha  known  polaa  could  ba 
polaa  of  tha  driving  function  which  ara  known  from  knowladga  of  tha  Laplaca 
tranaform  of  tha  driving  function,  ayatam  polaa  which  ara  known  from  pravloua 
Prony  analyala  or  othar  tachniquaa,  or  polaa  introduced  to  modal  tha  nolaa. 

It  la  wall  known  that  for  cartaln  data  aata  Prony' a  algorithm  haa 
aotua  difficulty  In  extracting  tha  trua  polaa  that  ara  contained  in  tha 
data.  Thla  problem  ia  generally  related  to  tha  nolaa  in  tha  data  but  will 
not  ba  dlacuaaad  hare.  Any  mathod  by  which  tha  accuracy  of  the  trua  polaa 
can  ba  incraaaad  will  ba  vary  uaaful. 

Logic  would  tall  ona  that  making  uaa  of  known  information  ahould 
incraaaa  tha  accuracy  of  a  calculation.  Hence,  if  uaa  la  mada  of  known 
polaa  ln  tha  data  it  would  aaam  raaaonable  to  aaauma  that  aoma  of  tha 
lnatabillty  ln  Prony' a  method  would  ba  alleviated.  Tha  proof  of  tha 
validity  of  thla  atatemant  raata  on  tha  actual  implementation  of  tha 
conatralnad  Prony  algorithm  and  tha  raaulta  compared  for  tha  aama  data 
analyzed  by  tha  unconatralnad  mathod. 

In  addition  to  aiding  ln  lncraaaing  tha  accuracy  of  tha  true  polaa 
it  might  ba  poaaibla  to  introduce  random  polaa  which  will  modal  tha  nolaa. 

In  that  mannar,  tha  polaa  that  ara  knoim  a  priori  ara  not  the  polaa  which 
wa  aaak.  Tha  difficulty  with  thla,  if  it  ahould  work,  ia  that  wo  need  to 
know  tha  rank  of  the  ayatam  ao  that  wo  can  introduce  the  proper  number  of 
nolaa  polaa. 


Another  asset  of  having  a  contrained  Prony's  method  is  that  it  can  be 
used  to  aid  in  deconvolution.  The  following  discussion  follows  directly 
from  work  performed  in  Reference  [A. 2]. 


Assume  that  a  system  is  excited  by  a  driving  function  which  can  be 
represented  in  the  form 


M 

G(t)  -  2 

J-l 


8j 


u(t). 


(H.l) 


It  io  presumed  here  that  the  driving  function  poles  s^  and  the  associated 
residues  g^  are  known.  These  can  be  determined  analytically  if  the  ana¬ 
lytical  form  of  the  driver  is  known  or  can  be  determined  from  a  Prony's 
method  fit  to  measurement  data  of  the  driving  function. 


Now  assume  that  the  response  of  the  system  to  the  driving  function 
G(t)  is 

M+N  s , t 

R(t)  -  2  r,  e  1  u(t)  (H.2) 

i-1  1 

and  that  the  impulse  response  of  the  system  is 
N  s  t 

H(t)  -  I  h.  e  k  .  (H.3) 

k-1  * 

Hence,  there  are  N  system  poles  and  residues  and  M  poles  in  the  driving 
function.  Of  course,  the  response  function  R(t)  can  be  written  as  the  con¬ 
volution  of  the  driving  function  G(t)  with  the  impulse  response  H(t).  That 
li» 

R(t)  -  H(tt)*G(t)  ,  ( H.4) 


where  *  denotee  convolution. 


What  is  usually  desired  is  to  obtain  tha  N  system  poles  and  residues 

of  (H.3).  This  must  be  done  by  deconvolution  of  the  driving  function 

(H.1)  from  (H.2).  If  the  r^  and  s^  are  known,  the  MfN  values  of  Sj,  contain 

the  M  values  s . .  The  g,  and  the  s.  are  known,  then  the  N  values  of  h  and 
J  J  j  k 


ak  can  be  determined  analytically. 


It  can  be  shown  that  the  response  function  can  be  written  in  terms  of 
the  driving  function  and  the  impulse  response  as 


M  N  h.  s.t  N  M  g.  s 

■  I  <8.  I  r^r)  • 1  ■  S  M  r^-> . 

J-l  3  k-l  s  k  k-1  k  J-l  aJ  *k 


(H  .5) 


Hence,  the  response  function  residues  can  be  defined  as 
N  hk 

r4  "  84  I  i  -  1  ,  M 

1  J  k-1  8j  8k 


(H.6a) 


M  8j 

N  I  ,  for  i  -  M+l  ,  N-Hi  . 

1  *  j-1  *j  “k 


(H.6b) 


From  (H.6b1  the  residues  of  the  impulse  response  are  written  as 


hk-|3: 

j-i  3j  ** 


k  -  1  ,  N 


i  -  M+k 


(H.  7) 


Thus  the  res'  * tes  or  amplitudes  of  the  impulse  response  can  be  obtained  from 
the  known  va.  as  r^,  g^,  s^,  s^.  Of  course,  since  the  M  values  of  are 
known  and  the  M+N  values  of  s^  are  presumed  to  contain  the  M  values  of  , 
then  the  N  values  of  s^  can  be  determined  by  inspection. 


Prony's  method,  however,  will  not  necessarily  give  the  true  values 
of  the  driving  function  poles  in  the  extracted  poles  of  the  response 
function.  Hence  a  constrained  Prony's  method  is  required  so  that  analyti¬ 
cal  deconvolution  of  (H.7)  cau  be  obtained. 

Work  performed  so  far  implements  Method  1  which  is  described  below. 

It  was  found  that  this  method  works  perfectly  as  long  as  no  noise  is 
present  in  the  signal.  As  soon  as  noise  is  introduced  in  the  signal  and 
a  least-squares  method  is  used  the  matrix  (H.17a)  is  corrupted  with  noise 
which  through  the  least-squares  process  also  corrupts  the  constrained 
parameter.  Experiments  have  confirmed  this  observation.  This  suggests 
that  the  constraining  Method  1  may  work  well  in  noisy  data  if  the  curve¬ 
fitting  Prony's  method  is  used.  This  has  not  been  tested  as  yet. 

Method  2,  forcing  the  polynomial  root  solver  to  find  certain  roots, 
has  also  not  been  tested.  This  is  because  if  the  coefficients  are  all 
corrupted  by  noise  through  the  least-squares  procedure  then  subtracting 
any  given  poles  out  of  the  polynomial  will  force  the  remaining  poles  to 
carry  the  burden  of  all  the  noise. 

H.2  Method  1 


In  the  implementation  of  Prony's  method,  an  Nth  order  polynomial  is 
solved  for  its  roots.  The  order  of  the  polynomial,  N,  is  the  number  of 
poles  being  sought  in  the  transient  data.  If  the  coefficients  of  the 
polynomial  are  denoted  as  a  then  the  polynomial  can  be  expressed  as  per 
Reference  H.l  as 

aQ  ♦  ajZ1  +  a2Z2  +  ...  +  a^ZN  -  0  (H.8) 
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js^aataastaaBBai 


where  is  usually  set  equal  to  unity.  The  N  roots,  Z^,  are  defined  as 


Zi  “  8 


(H.9) 


with  being  the  poles  sought,  and  At  being  the  time  step  size  used  in  the 
analysis. 


If  the  value  of  one  or  more  of  the  poles  is  known  -  that  is,  we  know 
some  of  the  -  then  the  Z^  can  be  substituted  into  (H.B).  For  example, 

if  s^  or  Z^  is  known,  then  (H.B)  can  be  written  as 


an  +  a.Z  +  a,Z  2  +  ...  +  a„Z  N  ■  0 
u  xx  xx  n  x 


(HIP) 


The  N+l  polynomial  coefficients  a  are  solved  for  in  Prony's  method  by 
solution  of  the  difference  equation 


■ _ 

^  “p  Xp+K  "  “XN+K’  ’  Y"M_N  * 


(HJ1) 


The  and  I^(R  are  the  samples  of  the  transient  signal  being  analyzed  and 
M  is  the  total  number  of  samples  being  used.  The  value  of  M  must  be  at  least 
equal  to  2N  to  give  N  sets  of  equations  in  the  N  unknowns  ctp.  However,  if 
the  value  of  a  pole  is  known,  then  one  of  the  N  aquations  can  be  equation 
H.10  and  N-l  equations  of  the  form  of  H.ll  can  be  written. 


If  L  poles  are  known  a  priori,  the  L  equations  of  the  form  of  (H.10) 


can  be  written  as 
N-l 

S  a  zj  ■ 

p-0  *  1 


-  Z  ,  l  ■  1 ,  2 ,  . . . ,  L 

& 


(H.12) 


and  y»M-N-L  equations  can  be  written  in  the  form  of  (H.ll) 

N-l 

p-o  °p  Ip+K  "  ~Ik+n  ’  k"0,;L . Y_1  ’  Y"M”N“L  •  (H.13) 

Hence  there  are  still  y-M-N  total  equations  to  solve  for  the  N  values  of  a 
however  the  system  is  constrained  by  the  knowledge  of  the  location  of  L  poles. 
As  is  usually  done  in  Prony's  method,  if  M-2N  then  the  set  of  equations  is 
inverted  and  solved.  If  M>2N  then  a  pseudo-inverse  procedure  is  used. 


Using  the  matrix  notation  of  Reference  H.l,  if  M-2N  and  L  poles  are  known, 
then  we  solve  the  equation 


AB  -  C 


where  A  is  a  square  matrix  defined  as 


A 


‘1 

,2 


•  •  •  Z 

•  •  •  Z 


N-l 

1 

N-l 
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.N-l 


...  I 


N-l 


...  I 
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(H.  14) 


(H.  15a) 
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and  B  and  C  are  vectors  defined  as 


“n-i 


(H. 15b) 


(H. 15b) 


L2N-1-L 


If  M>2N  and  L  poles  are  known,  then  the  solution  takes  on  a  pseudo-inverse 
or  least-squares  form  as 


T  T 
A  AB  ■  A  C 


(H.  16) 


Where  A  la  the  transpose  of  the  matrix  A  and  A  la  now  a  rectangular  matrix 
of  the  form 

i  z1  z 2  ...  zj’1 

l  z  z2  2n-i 

A  i>2  U  2  •••  Z2 


M-N-l-L 


...  I\ 


...  I. 


...  I 


M-L-2 


(H. 17a) 


(H. 17b) 


Consider  the  simple  example  of  a  two  pole  system 

st  st 

f(t)  -  2«  1  +  3s  L 


where  s1  ■  0  and  s2  •  -A.  Assume  that  s^  Is  known,  thus  giving  Z ^ 
Also,  let  f(t)  be  sampled  at  At  -  0.2  seconds,  giving  a  data  set  as 


2.6057 


2.2722 


2.1223 


Using  the  constrained  Prony's  method  in  the  square  system  form  of  Equation 
(H.14)  and  (H.15)  yields 


or 


1 

1 

ao 

"l 

5 

3.3480 

-°1- 

_2. 6057 

The  solution  of  this  sst  of  squations  givss  ■  1.4493  and  aQ  ■  0.4493  so 
that  tha  polynomial  can  ba  writtsn  as 

Z2  -  1.4493  Z  +  0.4493  -  0  . 

Tha  roots  of  this  polynomial  ara  Z^  ■  1.0  and  Z 2  •  0.4493,  giving  polas  of 
8^^  •  0.0  and  s2  ■  4.0003.  Tha  arror  in  ia  dua  to  truncation  arror.  Nota 
that  tha  constralnad  pola  was  ratumad  axactly. 

H.3  Mathod  2 

Anothar  approach  which  can  ba  usad  to  constrain  Prony's  mathod  to  usa 
information  about  known  polas  is  tha  modification  of  tha  polynomial  root 
finding  routina  MULLER.  Onca  tha  unconstrainad  Prony's  mathod  calculatas 
tha  coafficiants  of  tha  polynomial  (H.8)  and  if  soma  of  tha  roots  of  tha 
polynomial  ara  known  a  priori,  than  tha  locations  of  thosa  roots  can  ba 
passad  to  MULLER  and  it  will  not  hava  to  saarch  for  thosa  roots.  Sines  tha 
root  finding  routina  la  vary  tlma  consuming,  tha  knowlsdgs  of  tha  location 
of  any  of  tha  roots  will  presumably  save  computation  time. 

A  possible  flaw  with  this  approach  is  that  MULLER  is  forced  to  pre¬ 
sume  roots  of  the  polynomial  when  those  exact  roots  may  not  ba  contained  in 
the  polynomial.  That  is,  if  the  polynomial  has  not  been  constrained  to 
contain  the  known  roots,  as  per  Section  H.2,  then  the  known  roots  will 


not  necessarily  be  contained  In  It.  Forcing  an  unconstrained  polynomial 
to  have  certain  roots  could  grossly  perturb  the  location  of  the  other  roots 
being  sought. 

H.A  Method  3 

The  obvious  solution  to  the  flaw  presented  In  Section  H. 3  la  to  use 
both  the  methods  of  H.2  and  H.3  simultaneously.  That  1st  the  polynomial 
is  constrained  to  contain  the  known  poles  as  outlined  In  Section  H.2.  Then 
the  polynomial  root  finding  routine  MULLER  can  be  modified  to  extract  the 
known  roots  from  the  polynomial  before  It  begins  Its  search  for  the  unknown 
roots. 

It  Is  felt  that  this  approach  will  give  the  best  accuracy  and  will 
spaed  up  the  calculations  since  the  order  of  the  polynomial  Is  effectively 
reduced. 


:^i  '■n  ^:j:^,rw',^?^^uaiggp£r.,r,rv  ~ J 


Rafarancaa 


H.l  Van  Blaricum,  M.L.,  Tachnlguaa  for  Extracting  tha  Complax  Raaonancaa 
of  a  S vat am  Diractlv  from  lta  Tranaiant  Rasponaa.  Ph.D.  Diaaartation, 
Dacambar  1975 ,  Unlvaralcy  of  Illlnoia,  Dapartoant  of  E.E. 

H.2  A.  J.  Poggio,  M.L.  Van  Blaricum,  E.K.  Millar  and  R.  Mittra,  "Evaluation 
of  a  Procaaaing  Tachniqua  for  Tranaiant  Data,"  IEEE  Trana.  on  Ant,  and 
Prop. ,  pp.  165-173,  January  1973. 


H-ll 


Hi; 


tow  v  P™winn  j*"5"",*rTORTO,i^prtn rasiiRiirwrw^p^TTn  ^nif'njii  r»  n  i 
+**  '■«■*'»' . . . . *m**mJ*mm  ^-'■-  -'I 


ETOF^WW:  ,|PF  *’»i!*v»m»fcj  M»i ! 


APPENDIX  I 

REVERSING  THE  WAVEFORM  IN  TIME  TO  ELIMINATE 
EXTRANEOUS  RESONANCES 


Reversing  the  waveform  in  tiaa  depend*  on  Che  atatietical  properties 
of  Che  noise  which  do  not  change  when  Che  waveform  is  reversed.  The  poles 
chat  result  from  the  noise,  therefore,  are  not  altered  by  reversal  in  time 
while  the  true  poles  are  negated  or  flipped  through  the  origin. 

If  the  noise  level  is  high  the  noise  poles  only  approximately  remain 
the  same  under  time  reversal  and  the  true  poles  only  approximately  reflect 
through  the  origin.  If  time  reversal  is  attempted  for  curve-fitting  Prony's 
method  it  is  found  that  all  poles  flip  precisely.  Therefore,  this  method 
is  not  effective  for  curve-fitting  Prony's  method. 

(Curve-fitting  Prony's  method  uses  the  solution  to  the  Inhomogeneous  system 
QX*q,  in  the  notation  of  Volume  1,  where  Q  is  a  square  matrix.) 

The  waveform  of  Figure  1.1  was  used  in  a  numerical  example  of  the  timer" 
reversal  method.  Least-squares  Prony's  method  was  applied  to  the  waveform. 
Estimates  at  the  S-Plana  poles  were  found  using  the  inhomogeneous  solution 
which  is  defined  in  Volume  I,  Section  2  as 

\  -  [qpQf1  <fq. 

The  waveform  consists  of  100  samples  and  was  corrupted  with  uncorrelated, 
Gauaslan-distributed  noise  with  a  standard  deviation  of  0.1.  Figure  1.2 
displays  the  poles  obtained  from  the  forward  and  reversed  waveforms.  The 
dimensions  of  the  (NXn) -dimensional  matrix  Q  in  this  example  are  M-76 
and  r.**24.  The  estimates  of  the  true  poles  are  quite  accurate.  Figure  1.3 
displays  the  poles  obtained  for  the  same  waveform  but  different  matrix 
dimensions.  The  matrix  dimensions  are  M*60  and  n"40.  In  this  case,  the 
extraneous  poles  do  not  all  remain  in  the  left-half  plane  which  is 
attributed  to  using  an  overly  square  matrix. 
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If  the  matrix  Is  too  "thin",  of  M/n  »1,  the  estimates  of  the  true 
poles  become  Inaccurate.  If  the  matrix  is  too  "fat",  or  M/n  »  1,  the 
extraneous  poles  do  not  all  remain  in  the  left-half  plane.  The  matrix 
shape  where  M/n  «  .1  appears  to  be  the  best  shape  for  good  estimates  and 
keeping  the  extraneous  poles  in  the  left-half  plana. 
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Therefore,  It  appears  that  if  we  choose  M/n  *  3,  then  time 
reversal  can  be  used  to  distinguish  between  the  true  poles  and  the 
extraneous  poles  in  relatively  high  noise  levels. 

The  above  two  examples  show  that  there  is  a  trend  for  the 
noise  poles  in  the  least-squares  Prony's  method  to  occupy  a  higher 
frequency  portion  of  the  S-plane  and  to  be  highly  damped.  This  trend 
was  studied  by  letting  the  least-squares  Prony's  method  operate  on  a 
waveform  consisting  of  only  Gaussian-dlstributed  uncorrelated  noise  - 
no  signal.  The  poles  resulting  from  this  example  are  highly  damped  and 
evenly  distributed  in  the  z-plane  approximately  around  a  circle  within 
the  unit  circle. 


This  behavior  can  be  explained  by  the  fact  that  the  polynomial 

coefficients,  except  which  is  set  to  one,  all  tend  to  zero  for  the 

least-squares  method  operating  on  uncorrelated  noise.  The  polynomial 

N 

then  tends  toward  z  ■  e  where  e  tends  to  zero.  The  roots  of  this 
polynomial,  z^,  have  magnitudes 


It  then  follows  that  the  poles  should  be  highly  damped,  since 
|z^|  +  0  if  | c |  *0,  and  that  the  poles  are  evenly  distributed  about  the 
z-plane.  However,  for  curve-fitting  Prony's  method  the  a^,  i  *  0,  ..., 
N-l  do  not  ten  1  to  zero  and  indeed  this  phenomenon  is  not  observed  in 
tests  using  cui ve-fitting  Prony's  method  on  all  noise. 
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This  behavior  of  the  noise  poles  in  least-squares  Prony  seems  to 
remain  approximately  the  same  even  if  the  waveform  is  not  entirely  noise. 

The  trend  that  is  seen  is  that  the  noise  poles  occupy  the  higher  frequencies, 
are  highly  damped  and  evenly  distributed  between  the  higher  frequencies; 
while  the  true  poles  are  approximately  at  their  uncorrupted  locations 
and  occupy  the  lower  frequencies.  That  is,  it  appears  as  though  the  noise 
poles  are  "crowded"  away  from  the  lower  frequencies  or,  perhaps  more 
accurately,  the  lower  frequency  noise  poles  become  lower  frequency  signal 
poles  when  the  waveform  is  no  longer  entirely  noise. 


